Structural transition and re-emergence of iron's total electron spin in (Mg,Fe)O at ultrahigh pressure

Fe-bearing MgO [(Mg1−xFex)O] is considered a major constituent of terrestrial exoplanets. Crystallizing in the B1 structure in the Earth’s lower mantle, (Mg1−xFex)O undergoes a high-spin (S = 2) to low-spin (S = 0) transition at ∼45 GPa, accompanied by anomalous changes of this mineral’s physical properties, while the intermediate-spin (S = 1) state has not been observed. In this work, we investigate (Mg1−xFex)O (x ≤ 0.25) up to 1.8 TPa via first-principles calculations. Our calculations indicate that (Mg1−xFex)O undergoes a simultaneous structural and spin transition at ∼0.6 TPa, from the B1 phase low-spin state to the B2 phase intermediate-spin state, with Fe’s total electron spin S re-emerging from 0 to 1 at ultrahigh pressure. Upon further compression, an intermediate-to-low spin transition occurs in the B2 phase. Depending on the Fe concentration (x), metal–insulator transition and rhombohedral distortions can also occur in the B2 phase. These results suggest that Fe and spin transition may affect planetary interiors over a vast pressure range.

Authors, who are experts of this system, having studied already the structural/magnetic transition it undergoes in its B1 phase while compressed in the lower Earth's mantle pressure range, report a structural phase transition between B1 and B2 phases (accompanied by a transition from low to intermediate spin for Fe) and a second one in the B2 phase from IS to LS. The first of these is quite remarkable: in fact it makes the system re-acquire a spin (S = 1 of IS B2) from a state at S = 0 (LS B1), which goes against the tenet that pressure suppresses magnetization. In addition, for higher concentrations of Fe (x = 0.25) it also corresponds to an insulator-to-metal transition since the B2 phase is metallic.
Given the quality of the work and the originality of the results I think the paper deserves publication, even in the present form. The only comment I have actually concerns the metallic state of B2. As far as I understand authors judge this phase to be metallic based just on the partial occupation of t2g states (in the undistorted structure) and on the degeneracy of the multiplets resulting from the structural distortions. However, using these arguments is quite dangerous when working with DFT and band theory. In fact, there are many Mott insulators that, due to symmetry, appear to have a metallic band structure, although they are obviously non conductive. Have authors tried to further lower the symmetry of the system (e.g. eliminating all passible residual degeneracies) to check whether the metallic behavior would be robust against it? Although I strongly recommend the authors to take this point into consideration and to possibly update their paper accordingly I leave it as an optional modification and recommend the publication of the manuscript as is.

Reviewer #3 (Remarks to the Author):
The authors conducted first-principles calculations (LDA+U) to investigate the spin states and crystal structures of ferropericlase (Mg,Fe)O at TPa pressure range. They present two major results on two ferropericlase compositions (12.5% and 25% Fe; Fp15 and Fp25): (1). In 12.5% ferropericlase, B1 to rhombohedral B2 phase transition and the intermediate spin to lower spin transition in the rB2 phase; (2). in 25% ferropericlase, B1-B2 transition and the IS-LS spin transition. These results are surprising and warrant a publication in Nature Communications if they can be verified. Like the authors mentioned in the manuscript, the B1-B2 transition is expected to occur at high pressure in MgO-FeO system. However, the rhombohedral B2 phase is totally new. The reemergence of the spin state and the IS-LS transition in the B2 phase are unexpected and have important implications to our understanding of the interiors of exoplanets.
The authors should provide an explanation to why the Fp15 and Fp25 compositions have such a distinct structural transition behavior. The difference in the amount of iron in these two systems is rather small so one would not expect to see such a major difference. MgO-FeO system is a solid solution system so it is just a mystery to me that these two compositions display drastically different phase transition sequences. As shown in Fig 3, the lattice parameters in the rhombohedral B2 change drastically with pressure. Why is that? It will be useful if the authors construct a composition-pressure diagram to depict the phase transition sequence using literature data. Comparison to literature data in the MgO-FeO system will also help justify their calculations.

Reviewer's Comment
The predicted behavior of B2 phase of (Mg,Fe)O appears to be strongly sensitive to Fe concentration. The stability field of intermediate spin B2 phase for 25% iron case is much narrower than that for 12.5% iron case. Also, B2 phase is metallic for high iron concentration whereas it is an insulator for low concentration or pure MgO. As shown in Fig. 5, the enthalpy differences between different B2 phases are small, particularly, for intermediate spin state. It is important to assess how sensitive this energetics with respect to different computational factors. For instance, the values of U though obtained in a self-consistent way may not be uniquely defined. So the question is what happens if different U values are used. While LDA was used in this work, GGA is usually considered to be more appropriate for transition metals and iron-bearing systems. Previous studies (e.g., Holmstrom Stixrude, PRL, 2016;Ghosh and Karki, Sci. Rep., 2016) have shown that the predicted high spin to low spin transition of iron in (Mg,Fe)O differs significantly between GGA and LDA with/out the Hubbard term. More tests along this direction will be helpful.

Author's response
Before proceeding to further discussion, we would like to clarify that we had already performed extensive test calculations on (Mg1−xFex)O with various combinations of the functional (LDA/GGA) and Hubbard U parameter, as requested by the Reviewer. In our previous submission, GGA+U calculations with U = 6, 7, and 8 eV were already included in Supplementary Information (SI) but mislabeled as LDA+U. We apologize for such a mistake and the confusion it may have caused. In this revision, the GGA+U calculations ( Fig. S2) have been correctly labeled and also revised by including more data points for the Birch-Murnaghan EoS fitting. Newly added is Fig. S3, which shows additional LDA+U calculations with U = 10 and 14 eV. The procedure outlined in Figs. S2 and S3 is actually the common practice of DFT+U calculations: Treating U as a fitting parameter, performing calculations with various trial U and see how the calculation results are affected. Overall, the results obtained from the DFT+U test calculations with trial U (Figs. S2 and S3) are consistent with the LDA+Usc results shown in the main text (Fig. 6), as detailed below.
In Fig. S2, GGA+U calculations for (Mg1−xFex)O with trial U = 6, 7, and 8 eV are shown. For x = 0.125 [Figs. S2(a)−S2(c)], a transition from the B1 LS to the rB2 IS state occurs at ~0.63 TPa, followed by an IS-LS transition in the rB2 phase. For x = 0.25 [Figs. S2(d)-S2(f)], a transition from the B1 LS to the B2 IS state occurs at ~0.54 TPa, followed by an IS-LS transition in the B2 phase. Clearly, the predicted B1−(r)B2 structural-transition pressures are barely affected by U. In contrast, the predicted IS-LS spin-transition pressure in the (r)B2 phase significantly increases with U: As U increases from 6 to 8 eV, the IS−LS transition pressure increases from 1. The above DFT+U test calculations clearly indicate that the choice of the functional (LDA/GGA) and U parameter barely affects the predicted B1−(r)B2 transition pressure but significantly affects the IS−LS transition pressure in the (r)B2 phase. When U increases by 2 and 4 eV in GGA+U and LDA+U, respectively, the IS−LS transition pressure in the B2 phase increases by ~0.1 TPa (100 GPa) and > 0.2 TPa (200 GPa), respectively. Based on our own experience and the papers mentioned by the Reviewer, similar effects have been observed in other calculations for lower-mantle minerals: Increasing the trial U by 4 eV would increase the spin-transition pressure by several tens of GPa. Clearly, spin transition at high pressure cannot be accurately predicted using the trial U approach.
In principle, the Hubbard U parameter, which accounts for the on-site Coulomb interaction, should be computed from the first-principles, and the U value would depend on the iron valence/spin state and pressure. In this work, we compute the U parameter self-consistently (Usc) via a density functional perturbation theory (DFPT) approach (Ref. S6), which is equivalent to the linear response approach that we had used before (Refs. S3-S5). Within this approach, Usc can be determined unambiguously, opposed to the Reviewer's comment "the values of U though obtained in a selfconsistent way may not be uniquely defined". For years, we have extensively tested LDA+Usc and other DFT+U methods on Fe-bearing minerals of geophysical and/or geochemical importance, including B1 (Mg1−xFex)O (Ref. 22), Fe-bearing bridgmanite (Refs. 45,46) and post-perovskite MgSiO3 (Refs. 47,48), ferromagnesite (Mg1−xFex) CO3 (Refs. 49,51), and the new hexagonal aluminous (NAL) phase (Ref. 50). Throughout these works (Refs. 22,(45)(46)(47)(48)(49)(50)(51), we find that spin-transition pressure predicted by LDA+Usc is typically within 5-10 GPa around the experimental results, much more accurate than the trial U approach, which has a margin of error of up to tens or hundreds of GPa. The accuracy provided by LDA+Usc is of crucial importance for the predictive calculations on Fe-bearing minerals at ultrahigh pressure (TPa) relevant to exoplanet interior but challenging for experiments. We therefore continue adopting LDA+Usc in this work.
Comparing the trial DFT+U (Figs. S2 and S3) and LDA+Usc results (Fig. 6), the predicted B1−(r)B2 transition pressures are nearly the same, confirming the conclusion drawn from Figs. S2 and S3: B1−(r)B2 transition is barely affected by the choice of DFT functional and U. As to the IS−LS transition in the (r)B2 phase, the results of GGA+U with U = 7 [Figs. S2(b) and S2(e)] are in best agreement with the LDA+Usc results (Fig. 6). Figure S2 is correctly labeled as GGA+U and has also been revised by including more data points for Birch-Murnaghan EoS fitting. Figure S3 is newly added to show additional LDA+U results. Discussions for Figs. S2 and S3 (Sec. S3 in SI) have also been rewritten, similar to our response to the Reviewer.

Reviewer's Comment
The results analysis and discussion as presented in the paper look mostly fine. However, it is not clear what implications are, particularly, for planetary interiors. Can the authors be more specific in terms of density and composition.

Author's response
We are glad that the Reviewer considers the analysis and discussion as presented in the paper look mostly fine. Below we address the point raised by the Reviewer, about the implications for planetary interiors in terms of volume/density and composition.
(ii) Composition The interplay between the composition and spin transition can be complicated. The spin transition can be affected by the Fe concentration, and vice versa. For example, in the Earth's mantle condition, Fe partitioning between B1 (Mg1−xFex)O, Fe-bearing MgSiO3 bridgmanite, and postperovskite varies with pressure, temperature, and even the Fe valence/spin state (Refs. 52-54). Therefore, we can expect Fe concentration in B2 (Mg1−xFex)O varing with the depth or the interior region in exoplanet interiors, due to the variation of Fe partitioning between B2 (Mg1−xFex)O and other minerals phases, including post-perovskite and/or high-pressure silicates (Refs. 15,16). At the moment, however, the actual mineral phases and Fe partitioning in exoplanet interiors remain unknown. Therefore, we can at best discuss the possible effects of composition (Fe concentration) on spin transition as below, instead of the possible effects of spin transition on the composition.
By carefully comparing Figs. 6(a) and 6(b), we can notice that spin, structural, and metal-insulator transition of B2 (Mg1−xFex)O can be induced by the change of Fe concentration (x). As can be observed from

Change of manuscript
We have revised our manuscript based on the above discussion. The geophysical (or planetary) implications in terms of volume/density are included in Lines 200-211; the implications in terms of composition are included in Lines 167-185.
The study focusses specifically on two compositions of (M_1-x Fe_x)O: x = 0.125 and 0.250, sufficient to capture the effects of the growing concentration of Fe with increasing depth on the considered structural and elastic properties. For both concentrations the manuscript discusses the behavior of two structural phases, relevant in the considered pressure range: B1 and B2. The authors discuss thoroughly their electronic structure, symmetry and deformations, linking the structural response to pressure to the modifications the ordering of the Fe d states undergo when the crystal field changes. The effects of finite (high) temperatures are partially taken into account in the vicinity of phase transitions through a thermodynamic model based on the entropy of mixing of various phases at different spin (while the vibrational entropy is neglected in this work).
The manuscript is well structured and organized, written in a clear way, easy to read. The physics it discusses is sound and the results are original and quite important for the scientific community.
Authors, who are experts of this system, having studied already the structural/magnetic transition it undergoes in its B1 phase while compressed in the lower Earth's mantle pressure range, report a structural phase transition between B1 and B2 phases (accompanied by a transition from low to intermediate spin for Fe) and a second one in the B2 phase from IS to LS. The first of these is quite remarkable: in fact it makes the system re-acquire a spin (S = 1 of IS B2) from a state at S = 0 (LS B1), which goes against the tenet that pressure suppresses magnetization. In addition, for higher concentrations of Fe (x = 0.25) it also corresponds to an insulator-to-metal transition since the B2 phase is metallic.

Reviewer's Comment
Given the quality of the work and the originality of the results I think the paper deserves publication, even in the present form. The only comment I have actually concerns the metallic state of B2. As far as I understand authors judge this phase to be metallic based just on the partial occupation of t2g states (in the undistorted structure) and on the degeneracy of the multiplets resulting from the structural distortions. However, using these arguments is quite dangerous when working with DFT and band theory. In fact, there are many Mott insulators that, due to symmetry, appear to have a metallic band structure, although they are obviously non conductive. Have authors tried to further lower the symmetry of the system (e.g. eliminating all passible residual degeneracies) to check whether the metallic behavior would be robust against it?
Although I strongly recommend the authors to take this point into consideration and to possibly update their paper accordingly I leave it as an optional modification and recommend the publication of the manuscript as is.

Author's response
We thank the Reviewer for the strong recommendation "Given the quality of the work and the originality of the results I think the paper deserves publication, even in the present form". Throughout the Reviewer's report, only one point is raised: "The only comment I have actually concerns the metallic state of B2. ...Have authors tried to further lower the symmetry of the system (e.g. eliminating all passible residual degeneracies) to check whether the metallic behavior would be robust against it?". At the end of the report, however, the Reviewer decided to "leave it as an optional modification and recommend the publication of the manuscript as is".
To further improve our manuscript, despite that our response and modification for the manuscript are not required, we address the point raised by the Reviewer and revise our manuscript accordingly. It is well known that B1 (Mg1-xFex)O is insulating in the low Fe concentration regime (x £ 0.25), and the insulating high-spin (HS, S = 2) state is achieved via tetragonal Jahn-Teller (J-T) distortion, as further discussed in the next few paragraphs. Intuitively, we had also expected that metallic cubic B2 (Mg1-xFex)O with x £ 0.25 can be further stabilized and become insulating via distortion. Indeed, for x = 0.125, rhombohedrally-distorted rB2 (Mg0.875Fe0.125)O is insulating, as discussed in Lines 100−118 (Fig. 4) in the main text. Likewise, for cubic B2 (Mg0.75Fe0.25)O, we have also lowered its symmetry by applying rhombohedral compression/elongation for the IS/LS state, with the 3d orbital occupations guided accordingly. To our surprise, the entire crystal and the FeO8 polyhedra resumed cubic symmetry within a few steps of structural optimization. Consequently, the t2g orbitals are partially occupied, forming a partially-filled t2g band across the Fermi energy, and B2 (Mg0.75Fe0.25)O remains metallic. We have also tested tetragonal distortion, but after structural optimization, cubic symmetry and thus metallicity persisted. Effects of magnetic ordering were also investigated in these calculations. We find ferromagnetic (FM) state more energetically favorable than antiferromagnetic (AFM) state. In short, metallic cubic B2 (Mg0.75Fe0.25)O is very robust: No matter what we do, the cubic symmetry and metallicity persist. Descriptions of the aforementioned tests are included in the main text (Lines 125−127 and 60−62).
Next, we further clarify the relation between cubic symmetry and metallicity in the (Mg1-xFex)O system with x £ 0.25. We fully agree with the Reviewer that standard band theory and standard DFT methods , we broke the cubic symmetry by manually imposing tetragonal J-T distortion (stretching the Fe-O bonds on the xy plane), followed by structural optimization. After optimization, the structure remains tetragonal, and the energy is lower than the cubic symmetry. The three t2g orbitals split into to a singlet (dxy) and a doublet (dxz+ dyz), and the spin-down dxy orbital is occupied by one spin-down electron while the spin-down dxz and dzy orbitals are unoccupied [inset in Fig. R(b)]. Consequently, the spin-down t2g band splits into a completely filled dxy band below the Fermi level (0 eV) and a completely empty dxz+dyz band above 0 eV, as shown in Figs Fig. 4. For x = 0.25, however, the cubic symmetry is very robust (as described in the 2nd paragraph of our response). The main reason is that in B2 (Mg0.75Fe0.25)O, a 3D network of corner-sharing FeO8 cubes is formed, which suppress the rhombohedral distortion and further stabilizes the cubic structure, as detailed in Lines 119−141 (Fig. 5).

Change of manuscript
To better explain the metallicity of cubic B2 (Mg1-xFex)O, we have moved Fig. 3 from SI to the main text to discuss the metallicity of cubic B2 (Mg0.875Fe0.125)O, see Lines 79−99. We have also rewritten the discussion for B2 (Mg0.75Fe0.25)O (Fig. 5), see Lines 119−141. Our response to the Reviewer's comment about lowering the symmetry of cubic B2 (Mg0.75Fe0.25)O is included in the same part of the revision, see Lines 125−127.

Author's response
We thank the Reviewer for commenting our manuscript "surprising and warrant a publication in Nature Communications" and "have important implications to our understanding of the interiors of exoplanets". In the meantime, we would like to clarify one point: The work presented in this manuscript is a predictive calculation. To our knowledge, experiments for the MgO-FeO system at ultrahigh pressure have only been conducted on the end members MgO (Refs. 1-4 cited in the main text) and FeO (Refs. 43,44), not (yet) on Febearing (Mg1-xFex)O with 0 < x £ 0.25. Moreover, probing Fe spin states via currently available experimental techniques (e.g. x-ray emission spectroscopy and Mössbauer spectroscopy) at ultrahigh pressure has been highly challenging. Therefore, expecting or requesting for experimental verifications at this point would be unrealistic. Nevertheless, as described in Lines 49-57 of the main text, the computational method adopted in this work, LDA+Usc, has been shown to provide accurate prediction on spin transition of Fe-bearing minerals (Refs. 22,(42)(43)(44)(45)(46)(47)(48)(49)(50)(51). Moreover, to examine the robustness of our LDA+Usc results (Fig. 6), we have also performed GGA+U calculations (U = 6, 7, and 8 eV) and obtained nearly the same results, as detailed in Sec. S3 (see Fig. S2) of Supplementary  Information (SI). Therefore, the main conclusion of this paper, including the simultaneous LS-IS and B1-(r)B2 structural transition, and the IS-LS transition in the (r)B2 phase, should be reliable.

Reviewer's Comment
The authors should provide an explanation to why the Fp15 and Fp25 compositions have such a distinct structural transition behavior. The difference in the amount of iron in these two systems is rather small so one would not expect to see such a major difference. MgO-FeO system is a solid solution system so it is just a mystery to me that these two compositions display drastically different phase transition sequences. As shown in Fig 3, the lattice parameters in the rhombohedral B2 change drastically with pressure. Why is that?

Author's response
This comment from the Reviewer basically concerns the entire paper, namely, the unique properties of B2 (Mg1-xFex)O, which include: (1) The structural and electronic properties are sensitive to the Fe concentration, even at low concentration (x ≤ 0.25), see Figs. 3−5 in the main text, and (2) the effects of the composition on the structural and spin transitions of (Mg1-xFex)O, see Fig. 6. Prior to this paper, we had also been expecting the properties of B2 (Mg1-xFex)O to remain similar throughout 0 < x ≤ 0.25. To our surprise, properties of B2 (Mg1-xFex)O changes drastically as the Fe concentration increases from x = 0.125 (referred to as "Fp15" by the Reviewer, although it should be "Fp12") to x = 0.25 (referred to as "Fp25" by the Reviewer). Given this unexpected result, the main focus of this paper is therefore to analyze, explain, and discuss the unique properties of the B2 phase in detail. Also, despite the aforementioned composition dependence, both compositions exhibit resemblance in the sequence of transitions: B1 LS → rB2 IS → rB2 LS for x = 0.125, and B1 LS → B2 IS → B2 LS for x = 0.25 (Fig. 6). For some reason, however, these major points in the original manuscript were not successfully conveyed to the Reviewer. Here, we further clarify these points and also revise our manuscript for better clarity, as described below.
Before discussing the newly reported (r)B2 (Mg1-xFex)O in this paper, we first use the well-known B1 phase as an example to explain the relation between lattice distortion and insulating state in the (Mg1-xFex)O system with x ≤ 0.25. In the B1 phase, Fe substitute Mg in the 6-coordinate site, forming FeO6 octahedra, with the t2g orbitals having lower energy than the eg orbitals, as shown in Fig. S1 of Supplementary Information (SI). At low pressure, B1 (Mg1-xFex)O with x £ 0.25 is known to be an insulator with high-spin (HS, S = 2) Fe 2+ . This HS insulating state is stabilized via tetragonal Jahn-Teller (J-T) distortion, with longer Fe-O bonds on the xy plane and slightly longer lattice parameters in the x and y direction. Without J-T distortion, the HS state would be metallic and has higher energy. After structural optimization, the structure remains tetragonal, and the energy is lower than the cubic symmetry. With J-T distortion, the three t2g orbitals split into a singlet (dxy) with lower energy and a doublet (dxz and dyz) with higher energy [inset in Fig. R(b)]. For the HS state, the spin-down dxy orbital is occupied by one spin-down electron, and the spin-down dxz and dzy orbitals are unoccupied [inset in Fig. R(b)]. Consequently, the spin-down t2g band splits into a completely filled dxy band and a completely empty dxz+dyz band, as shown in Figs. R(b) and R(d). Evidently, for B1 (Mg1-xFex)O with x £ 0.25, the incorrect metallic HS state is a consequence of cubic symmetry, and the correct insulating HS state is achieved via distortion.
Similarly, cubic B2 (Mg1-xFex)O with x ≤ 0.25 is also metallic, with a partially-filled t2g band. To achieve an insulating state in this phase, lattice distortion would be necessary. Indeed, for x = 0.125, we find cubic B2 (Mg0.875Fe0.125)O dynamically unstable, and the IS/LS state is stabilized via rhombohedral compression/elongation, referred to as rB2 IS/LS (Mg0.875Fe0.125)O, as detailed in Lines 79-99 (Fig. 3) and 100-118 (Fig. 4) in the main text. (Note: Figure 3 was formerly placed in SI of the original manuscript.) For x = 0.25, to our surprise, the cubic structure is very robust: When we applied rhombohedral distortion on B2 (Mg0.75Fe0.25)O, the crystal structure resumed cubic symmetry within a few steps of the structural optimization. Furthermore, cubic B2 (Mg0.75Fe0.25)O is dynamically stable.
The main reason is that in B2 (Mg0.75Fe0.25)O, a three-dimensional (3D) network of corner-sharing FeO8 is formed, which suppresses the rhombohedral distortion and further stabilizes the cubic structure, as detailed in Lines 119-141 (Fig. 5). Within the cubic structure, B2 (Mg0.75Fe0.25)O is metallic, regardless the spin state. For the Reviewer's convenience, these results are further described below (but briefer than the manuscript  Fig. 4]. (Note: Figure 4 in the revised manuscript was formerly Fig. 3 in the original manuscript.) Briefly speaking, the IS/LS state is stabilized via rhombohedral compression/elongation, with shortened/stretched Fe-O bonds along the [111] direction, and the rhombohedral angle (a) greater/smaller than 90° [Figs. 4(a) and 4(f)]. In the crystal field with rhombohedral compression/elongation, the t2g orbitals split into a singlet a1g and a doublet e'g, and the a1g orbital has higher/lower energy than the e'g orbitals, due to the shortened/stretched Fe-O bonds along the [111] direction [Figs. 4(c) and 4(h)]. For the IS state, the spin-up e'g orbitals are completely occupied by two spin-up electrons, and the spin-up a1g orbital is completely unoccupied [Fig. 4(c)]. For the LS state, the a1g orbital is completely occupied, and the e'g orbitals are completely unoccupied [ Fig. 4(h)]. Consequently, a gap is opened between the a1g and e'g bands [Figs. 4(e) and 4(j)]. Clearly, the insulating state of rB2 (Mg0.875Fe0.125)O is a direct consequence of rhombohedral distortion, similar to the insulating HS state in the B1 phase achieved via tetragonal J-T distortion (Fig. R).
To our surprise, in B2 (Mg0.75Fe0.25)O (namely, x = 0.25), rhombohedral distortion does not occur, and B2 (Mg0.75Fe0.25)O remains cubic (see Fig. 5). Briefly speaking, with x = 0.25, the Fe concentration is high enough to form a 3D network of corner-sharing FeO8 cubes; in contrast, for x ≤ 0.125, the FeO8 polyhedra are isolated (unconnected), as shown in Fig. 1 (bottom panels). For isolated FeO8 polyhedra, rhombohedral distortion is allowed and favored, as observed in rB2 (Mg0.875Fe0.125)O (Fig. 4). The connectivity of the 3D FeO8 network in B2 (Mg0.75Fe0.25)O, however, suppresses the rhombohedral distortion and further stabilizes the cubic structure. When performing structural optimization, we did apply rhombohedral distortion on B2 (Mg0.75Fe0.25)O, but within a few steps, the crystal structure and FeO8 polyhedra resumed cubic symmetry. Based on these tests, dynamical stability of cubic B2 (Mg0.75Fe0.25)O can thus be expected [ Fig. 5 Finally, the Reviewer seemed to misunderstand Fig. 4 (formerly Fig. 3 in the original manuscript). In Fig. 4, we plot the atomic structure of rB2 (Mg0.875Fe0.125)O at volume 55.310 Å 3 /f.u. (~1.07 TPa), simply to point out that the IS and LS states are stabilized via rhombohedral compression and elongation, respectively, and therefore, have different structural properties, as explained in the 4th paragraph of our response and also detailed in Lines 100-118. Figure 4 has nothing to do with spin transition. Spin transition is determined by comparing the enthalpy of these two states, as shown in Fig. 6(a): A crossing occurs at 1.438 TPa, indicating a transition from rB2 IS to rB2 LS state. Given that rB2 IS and LS states have slightly different structures, this transition is in fact a simultaneous spin and structural transition: rhombohedrally-compressed IS state ( > 90°) → rhombohedrally-elongated LS state ( < 90°). It is well known that spin transition is often accompanied by a change of structure, as different spin states have different lattice distortions (due to different orbital occupations).

Reviewer's Comment
It will be useful if the authors construct a composition-pressure diagram to depict the phase transition sequence using literature data. Comparison to literature data in the MgO-FeO system will also help justify their calculations.

Author's response
As mentioned in our response to the Reviewer's first comment, to our best knowledge, there is no experimental data in the literature to make a composition-pressure diagram for the MgO-FeO system at ultrahigh pressure. Currently, ultrahigh-pressure experiments for the MgO-FeO system have only been conducted on  and FeO (Refs. 43,44). As described in the main text, experiments have shown that the B1-B2 transition of MgO and FeO occurs at ~0.5 TPa (see Line 23) and ~0.25 TPa (see Line 44), respectively. In our calculation, we find the B1-B2 transition of MgO occurs at 0.535 TPa, consistent with previous calculations . As the Fe concentration increases from x = 0 (MgO) to x = 0.125, the B1-B2 transition pressure increases to 0.642 TPa [ Fig.  6(a)]. Remarkably, when the Fe concentration further increases to x = 0.25, the B1-B2 transition pressure decreases to 0.539 TPa [ Fig. 6(b)], indicating a trend that the B1-B2 transition pressure decreases with Fe concentration for x > 0.125. This result is consistent with experiments: a much lower B1-B2 transition pressure in FeO (~0.25 TPa) than in MgO (~0.5 TPa).

Change of manuscript
The above discussion is now included in the manuscript, see Lines 157-163 in the main text.